% Generate specifications for the 9 contracts presented in the experiment.  
    % --> stored in file 'ContractSpecs.mat'

Cost_grid = [0.002, 0.003, 0.005, 0.007];   C_num = length(Cost_grid);

%% financial market data

r = 0.03;           % constant annual interest rate, compounded continuously.
lambda = 0.20;      % likelihood of jump.
mu_d_RN = 0.0735;   % mean of diffusion process (under risk-neutral measure).
sigma_d = 0.14;     % volatility of diffusion process.
mu_j = -0.25;       % mean of each jump.
sigma_j = 0.10;     % volatility of magnitude of each jump.

%% Contract specs

num_con = 9;                    % number of contracts used in experiment.

con_type_vec(1:3) = "TDF";
con_type_vec(4:6) = "Buffer";
con_type_vec(7:9) = "Floor";

con_spec_vec = zeros(50,num_con);    % contains 'alpha' for TDF, 'B' for buffer-RILA and 'F' for floor-RILA.

% TDFs:
con_spec_vec(1:20,1) = 1;       con_spec_vec(21:50,1) = 0.95 - 0.01*(1:30);
con_spec_vec(:,2) = con_spec_vec(:,1) - 0.10;
con_spec_vec(:,3) = con_spec_vec(:,1) - 0.20;

% Buffer RILAs:
con_spec_vec(1:20,4) = 0;       con_spec_vec(21:50,4) =1/300 * (1:30);
con_spec_vec(:,5) = con_spec_vec(:,4) + 0.025;
con_spec_vec(:,6) = con_spec_vec(:,4) + 0.050;

% Floor RILAs:
con_spec_vec(1:20,7) = 0.25;  	con_spec_vec(21:50,7) = 0.25 - 1.9/300*(1:30);
con_spec_vec(:,8) = con_spec_vec(:,7) - 0.035;
con_spec_vec(:,9) = con_spec_vec(:,7) - 0.055;


%% Find Cap Rates of RILA contracts:

con_cap_vec = zeros(50,C_num,num_con);      % contains cap rates for RILA contracts.

% simulate log returns of index for 1 year:
N_sim = 1e+6;   % number of simulated paths.
num_jumps = poissrnd(lambda,N_sim,1);
logR = (mu_d_RN - sigma_d^2/2)*1 + num_jumps*mu_j + sqrt( sigma_d^2 * 1 + num_jumps * sigma_j^2 ) .* randn(N_sim,1);   % vector of annual log returns, following MJD model.

options = optimset('Display','none');

for t = 50:-1:1
    for c = 1:C_num
        c_bar_RILA = Cost_grid(c);
        for j = 1:num_con
            con_type = con_type_vec(j);
            if con_type == "TDF"
                % do nothing
            elseif con_type == "Buffer"
                con_cap_vec(t,c,j) = fsolve(@(C) ExpLoss("Buffer",con_spec_vec(t,j),C,r,1,0,c_bar_RILA,logR) , 0.1 , options);
            elseif con_type == "Floor"
                con_cap_vec(t,c,j) = fsolve(@(C) ExpLoss("Floor",con_spec_vec(t,j),C,r,1,0,c_bar_RILA,logR) , 0.1 , options);
            else
                disp('Contract type unspecified.')
            end
        end
    end
end

save('ContractSpecs.mat','Cost_grid','con_type_vec','con_spec_vec','con_cap_vec')


%% Function to help find fair RILA cap rates

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function EL = ExpLoss(D_type,D_level,Cap,r,tau,fee,c_bar,logR)

R_S = exp(logR) - 1;                                                % simulated annual net returns of market index.
A_t = 1000;
if D_type == "Buffer"
    B = D_level;
    R_A = min(R_S+B , 0) .* (R_S<0) + min(R_S,Cap) .* (R_S>0);      % vector (N_sim-by-1) of credited return to RILA account.
elseif D_type == "Floor"
    F = D_level;
    R_A = max(R_S , -F) .* (R_S<0) + min(R_S,Cap) .* (R_S>0);       % vector (N_sim-by-1) of credited return to RILA account.
else
    disp('Error. Downside protection type unspecified.')
    asdfadsf    % to produce error
end

A_tplustau = A_t * exp(-fee*tau) * (1+R_A);
L = A_tplustau * exp(-r*tau) + A_t * c_bar * tau - A_t;
EL = mean(L);
return
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%